function oz_snapshot,isnap,itype,var,indir=indir,dir_prefix=dir_prefix,file_prefix=file_prefix,ics=ics,verbose=verbose

;+
; NAME:
;       OZ_SNAPSHOT
;
; PURPOSE:
;       New routine to read in single variables for single particle
;       types from gadget hdf5 snapshots. Levers new hdf5 library
;       commands to circumvent invoking the awkward k-shell scripts
;       used by the 'snap.pro' Tom and I developed
;
; CATEGORY:
;       
;
; CALLING SEQUENCE:
;      snapshot(isnap,itype,var)
;
; INPUTS:
;       isnap - integer index of the snapshot desired
;       itype - integer index of the particle type required [0:5]
;       var   - string containing the name of the desired variable
;
; OPTIONAL INPUTS:
;       indir      - changes the default directory path from './'
;       dir_prefix - changes the default directory prefix from
;                    'snapshot_'
;       file_prefix - changes the default file prefix from 'snap_'
;       verbose     - gives more output info than usual
;       ics         - flags that snapshot is an IC file, so has no
;                     snapshot number, meaning isnap is ignored
;
; OUTPUTS:
;       structure containing the relevant data
;
; MODIFICATION HISTORY:
;       Started by Rob Crain, Swinburne, 22-Nov-2010
;-

; Default values
if keyword_set(indir)then begin
    my_dir = indir
endif else begin
    my_dir         = '.'
endelse
if keyword_set(dir_prefix)then begin
    my_dir_prefix = dir_prefix
endif else begin
    my_dir_prefix  = 'snapshot_'
endelse
if keyword_set(file_prefix) then begin
    my_file_prefix = file_prefix
endif else begin
    my_file_prefix = 'snap_'
endelse

; Warnings and errors
if keyword_set(ics) && keyword_set(verbose) then print,' [snapshot] flagged as ICs, so ignoring the input snapshot number'
if (n_params() LT 3) then goto, error
if (itype LT 0 || itype GT 5) then goto, error

; Write numerics into string
this_snap = string(isnap,format='(i3.3)')

; Try opening single file first
nfiles = 1
if keyword_set(ics)then begin
   infile = my_dir+'/'+my_dir_prefix+'/'+my_file_prefix+'.hdf5'
endif else infile = my_dir+'/'+my_dir_prefix+this_snap+'/'+my_file_prefix+this_snap+'.hdf5'
file_exists = file_test(infile)

; If file not found, first try the multifile option
if(file_exists EQ 0)then begin 
   if keyword_set(ics)then begin
      infile = my_dir+'/'+my_dir_prefix+'/'+my_file_prefix+'.0.hdf5'
   endif else infile = my_dir+'/'+my_dir_prefix+this_snap+'/'+my_file_prefix+this_snap+'.0.hdf5'
   file_exists = file_test(infile)

   ; If still not found, return with a -1 value
   if(file_exists EQ 0)then begin 
      if keyword_set(verbose)then print,'  [snapshot] cannot find snapshot in single or multi-file format' 
      if keyword_set(verbose)then print,'  [snapshot] ',infile
      return,-1
   endif 
endif

; Determine number of files in this multi-file snapshot
NumFiles = h5ra(infile,'Header','NumFilesPerSnapshot')
if keyword_set(verbose)then print,' [snapshot] NumFiles = '+strtrim(string(NumFiles),2)

; Determine Redshift
Redshift = h5ra(infile,'Header','Redshift')
if keyword_set(verbose)then print,' [snapshot] Redshift = '+strtrim(string(Redshift),2)

; Determine number of particles of each species
NumPart_Total = h5ra(infile,'Header','NumPart_Total')
if keyword_set(verbose)then $
   for jtype=0,5 do print,' [snapshot] NumPart_Total['+strtrim(string(jtype),2)+'] = '+strtrim(string(NumPart_Total[jtype]),2)


; loop over files to collect data
varname     = 'PartType'+string(itype,format='(i1.1)')+'/'+var
rank        = 0
ncur        = 0
running_tot = 0 
itrip       = 0
result      = -1
for ifile=0,NumFiles-1 do begin
   
   if(NumFiles eq 1) then begin
      if keyword_set(ics)then begin
         infile = my_dir+'/'+my_dir_prefix+'/'+my_file_prefix+'.hdf5'
      endif else infile = my_dir+'/'+my_dir_prefix+this_snap+'/'+my_file_prefix+this_snap+'.hdf5'
   endif else begin
      if keyword_set(ics)then begin
         infile = my_dir+'/'+my_dir_prefix+'/'+my_file_prefix+'.'+strtrim(string(ifile),2)+'.hdf5'
      endif else infile = my_dir+'/'+my_dir_prefix+this_snap+'/'+my_file_prefix+this_snap+'.'+strtrim(string(ifile),2)+'.hdf5'
   endelse

   ; Check for the file
   file_exists = file_test(infile)
   if(file_exists EQ 0) then begin
      if keyword_set(verbose)then print,' [snapshot] cannot find file',infile
      return,-1
   endif else begin
      tmp = h5rd(infile,varname)
      nn = size(tmp,/dimen)
      rank = n_elements(nn)
      if(rank gt 1)then begin
         if(itrip eq 0) then begin
                                ; use itrip here instead of if(ifile eq 0) because
                                ; sometimes the first file contains no particles
                                ; of a given type
            result = fltarr(nn[0],NumPart_Total[itype])
            itrip = 1
         endif
         npart  = nn[1]
      endif else begin
         if(itrip eq 0) then begin
            if(strcmp(var,'Coordinates') eq 1)then begin
                                ; this is pathological case of there
                                ; being just one particle in a file
                                ; it *has* happened, esp. for PartType5    
               result = fltarr(3,NumPart_Total[itype])
            endif else begin
               case var of
                  'ParticleIDs':           result = lonarr(NumPart_Total[itype])
                  'Lx_MEKAL_0p5_2p0':      result = dblarr(NumPart_Total[itype])
                  'Lx_APEC_0p5_2p0':       result = dblarr(NumPart_Total[itype])
                  'Lx_APEC_bol':           result = dblarr(NumPart_Total[itype])
                  'APEC/Lx_MEKAL_0p5_2p0': result = dblarr(NumPart_Total[itype])
                  'APEC/Lx_APEC_0p5_2p0':  result = dblarr(NumPart_Total[itype])
                  'APEC/Lx_APEC_bol':      result = dblarr(NumPart_Total[itype])
                  else:                    result = fltarr(NumPart_Total[itype])
               endcase
            endelse
            itrip = 1
         endif
         if(strcmp(var,'Coordinates') eq 1)then begin 
            npart  = nn[0]/3
            nn = lonarr(2)
            nn[0] = 3
            nn[1] = 1
            rank = 2
         endif else begin
            npart  = nn[0]
         endelse
         
      endelse

      running_tot = running_tot + npart
      if keyword_set(verbose) then print,' [snapshot] Read snapshot '+strtrim(string(infile),2)+ $
                                         ', contains '+strtrim(string(npart),2)+' particles.'
      if(rank gt 1)then begin
         for ic=0,nn[0]-1 do begin
            result[ic,ncur:ncur+npart-1] = tmp[ic,*]
         endfor
      endif else begin
         result[ncur:ncur+npart-1] = tmp[0:npart-1]	
      endelse
      ncur = ncur + npart
   endelse
endfor

return,result

error: print,' usage: function snapshot(isnap,PartType,Var,indir=indir,dir_prefix=dir_prefix,file_prefix=file_prefix)'
       print,' read variable Var for species PartType [0-5] from snapshot number isnap'
       print,' indir = indir : optionally use other than default directory ',my_dir
       print,' dir_prefix = dir_prefix : optionally use other than default snapshot',my_dir_prefix
       print,' file_prefix = file_prefix : optionally use other than default snapshot file ',my_file_prefix
       print,' option verbose'
return,-1

end
